#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          //////////////////

          char kappa_sea[1024];
          int conf_ini, jump;
          int  Nex;
          int Nmass,Nr;
	  int NOa2;
	  int Nf;
          double beta, Uplaq;
          char dir_input[1024];
          char dir_output[1024];
          char gauge_fixing[1024];
          char gauge_action[1024];



          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
            if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
            if(strncmp(argv[i],"-jump",5)==0)  sscanf(argv[i],"-jump=%d",&jump);
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);
            if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-Nr",3)==0)  sscanf(argv[i],"-Nr=%d",&Nr);
            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);

          }
          printf("==============================\n");
          printf("beta=%0.3f Uplaq=%0.5f\n",beta,Uplaq);
          printf("==============================");

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);




	// IO

        char file[1024];
	char Oa2print[1024];
	char z_output[1024];

       //////////////////////////////////
        // Read Zq and LAMBDA
        /////////////////////////////
        printf("====================\n");
        printf(" READ Zq and LAMBDA    \n");

        double ****Zq=malloc(sizeof(double***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zq[iclust][imass]=malloc(sizeof(double*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zq[iclust][imass][imom]=malloc(sizeof(double)*(NOa2));

	double *****Lambda=malloc(sizeof(double****)*(5));
        for (int igamma=0; igamma<5; igamma++) Lambda[igamma]=malloc(sizeof(double***)*(Nclust+1));
        for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust]=malloc(sizeof(double**)*(Nmass*Nmass));
        for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) Lambda[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom));
	for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Lambda[igamma][iclust][imass][imom]=malloc(sizeof(double)*(NOa2));

      double **aux_allmomenta=malloc(sizeof(double*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

    // Zq
    for (int imass=0; imass<Nmass; imass++)
    {
         sprintf(file,"Zq_m%d",imass);
         read_input_momentabymomenta(dir_output,file,aux_allmomenta);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imass][imom][0]=aux_allmomenta[imom][iclust];
    }

    // ZO 
   char gammaprint[8];
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
        {
         int imass=im2*Nmass+im1;
         sprintf(file,"Lambda%s_m%d%d",gammaprint,im1,im2);
         read_input_momentabymomenta(dir_output,file,aux_allmomenta);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust][imass][imom][0]=aux_allmomenta[imom][iclust];
       }
    }

     printf("====================================================\n");

     ///////////////////////////////
     // Compute O(a2g2) corrections
     ///////////////////////////////

     printf("====================================================\n");
     printf("Zq Z_S Z_V Z_P Z_T Z_A with O(a2g2) corrections \n ");


     //////////////// substracting terms     
     double ***CorOa2_Zq=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Zq[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Zq[iOa2][imom]=malloc(sizeof(double)*(1));
     double ***CorOa2_Lambda=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Lambda[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Lambda[iOa2][imom]=malloc(sizeof(double)*(5));

     //iOa2=0 -> uncorrected
     for (int imom=0; imom<Nmom; imom++){
     CorOa2_Zq[0][imom][0]=0.;
     for (int igamma=0; igamma<5 ; igamma++) CorOa2_Lambda[0][imom][igamma]=0.;
     }
     //iOa2=1 corrected
     //iOa2=2 corrected boosted
     double g2;
     for (int iOa2=1; iOa2<NOa2; iOa2++)
     {
        for (int imom=0; imom<Nmom; imom++)
        {
          if(iOa2==1) g2=6.0/beta;
          if(iOa2==2) g2=(6.0/beta)/Uplaq;
          correction_Oa2g2(CorOa2_Zq[iOa2][imom],CorOa2_Lambda[iOa2][imom],g2,SinP2[imom],SinP4[imom],nn[imom],gauge_action,gauge_fixing);
        }

     }
      ////////////////// Zq corrected
      for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++)
      for (int imass=0; imass<Nmass; imass++)
      {
         for (int iOa2=1; iOa2<NOa2; iOa2++) Zq[iclust][imass][imom][iOa2]=Zq[iclust][imass][imom][0]-CorOa2_Zq[iOa2][imom][0];
      }

      ///////////////// Lambda corrected
    for (int igamma=0; igamma<5 ; igamma++) 
    for (int iOa2=1; iOa2<NOa2; iOa2++)
    for (int imass=0; imass<Nmass*Nmass; imass++) 
    	for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++)
	  Lambda[igamma][iclust][imass][imom][iOa2]=Lambda[igamma][iclust][imass][imom][0]-CorOa2_Lambda[iOa2][imom][igamma];


      //////////////// ZS ZV ZP ZT ZA corrected
    double *****ZO=malloc(sizeof(double****)*(5));
    for (int igamma=0; igamma<5; igamma++) ZO[igamma]=malloc(sizeof(double***)*(Nclust+1));
    for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust]=malloc(sizeof(double**)*(Nmass*Nmass));
    for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) ZO[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom));
    for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) ZO[igamma][iclust][imass][imom]=malloc(sizeof(double)*(NOa2));


    for (int iOa2=0; iOa2<NOa2; iOa2++)
    for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++) 
    for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
    {
	int imass=im2*Nmass+im1;
	ZO[igamma][iclust][imass][imom][iOa2]=sqrt(Zq[iclust][im1][imom][iOa2]*Zq[iclust][im2][imom][iOa2])/Lambda[igamma][iclust][imass][imom][iOa2];
//	ZO[igamma][iclust][imass][imom][iOa2]=sqrt(Zq[iclust][im1][imom][0]*Zq[iclust][im2][imom][0])/Lambda[igamma][iclust][imass][imom][iOa2];
    }

    
    printf("====================================================\n\n\n");

        /////////////////
        // Write output
        ////////////////



     double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));


      printf(" WRITING OUTPUT IN %s\n",dir_output);

    ////////////////
    // Zq R-PROM
    ////////////////////
   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
    whichOa2(iOa2,Oa2print);
    for (int imass=0; imass<Nmass; imass++)
     {
         sprintf(z_output,"%s/Z_ALLP/Zq_%s_m%d.out",dir_output,Oa2print,imass);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq[iclust][imass][imom][iOa2];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2);

         sprintf(file,"Zq_%s_m%d",Oa2print,imass);
         write_output_momentabymomenta_real(dir_output,file,aux_allmomenta_real,SinP2);
     }
   }
   ///////////////////////////////
   /// ZO R-PROM & GAMMA-PROM
   ///////////////////////////////////
   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
   whichOa2(iOa2,Oa2print);
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
         {
         int imass=im2*Nmass+im1;

         sprintf(z_output,"%s/Z_ALLP/Z%s_%s_m%d%d.out",dir_output,gammaprint,Oa2print,im1,im2);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=ZO[igamma][iclust][imass][imom][iOa2];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2);

         sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
         write_output_momentabymomenta_real(dir_output,file,aux_allmomenta_real,SinP2);


          }
   }
   }

	free(Zq); free(Lambda);
    return 0;

}

